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Abstract 

We present a simulation algorithm that accurately propagates a 
molecule pair using large time steps without the need to invoke the 
full exact analytical solutions of the Smoluchowski diffusion equation. 
Because the proposed method only uses uniform and Gaussian random 
numbers, it allows for position updates that are two to three orders of 
magnitude faster than those of a corresponding scheme based on full 
solutions, while mantaining the same degree of accuracy. Neither sim¬ 
plifying nor ad hoc assumptions that are foreign to the underlying 
Smoluchowski theory are employed, instead, the algorithm faithfully 
incorporates the individual elements of the theoretical model. The 
method is flexible and applicable in 1, 2 and 3 dimensions, suggesting 
that it may find broad usage in various stochastic simulation algo¬ 
rithms. We demonstrate the algorithm for the case of a non-reactive, 
irreversible and reversible reacting molecule pair. 


1 Introduction 

Naive Brownian dynamics (BD) simulations notoriously suffer on only ineffi- 
cently modeling a molecule’s diffusive motion near boundaries. In fact, tiny 
time steps are in general necessary to resolve the diffusive behavior close 
to boundaries with acceptable accuracy, rendering the naive application of 
the BD scheme highly incapable of accounting for reaction-diffusion systems 
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over biologically relevant time scales. Many particle-based stochastic simu¬ 
lation algorithms seek remedy by employing exact analytical solutions of the 
Smoluchowski diffusion equation that incorporate suitable boundary condi¬ 
tions (BC) [gang EE Hamm El. Because the analytical solutions take 
into account the correct boundary behavior, much larger time steps can be 
used without sacrifying accuracy. 

Schematically, a particle-based algorithm describing bimolecular reac¬ 
tions consists of three major parts. First, it has to decide if an encounter 
took place during a simulation time step. Secondly, if there was an en¬ 
counter, did also a reaction occur? Thirdly, if no reaction occurred, the en¬ 
counter pair is propagated according to a full solution of the Smoluchowski 
equation that takes into account the BC. In this general scheme, the step 
that involves the propagation is the most time consuming one, it involves 
sampling from a probability density function (PDF) that is given by a se¬ 
ries expansion where the individual terms are represented by complicated, 
numerically difficult to evaluate, integrals. 

This is where the proposed method comes into play: It eliminates the 
need to sample from a complicated series expansion, without giving up on 
accuracy or large size of the time step. Importantly, only the Gaussian PDF 
has to be sampled from, which results in a tremendously faster execution 
of the propagation move. While other algorithms using Gaussian random 
numbers to update the particle positions have been suggested mm , our 
approach does not require rescaling of reaction rates based on macroscopic 
parameters (2) and is particularly close to the underlying physics as well as, 
at the same time, easy to implement. Because the position moves provide 
en passant the first-passage (FP) times, it facilitates the use of small-time 
approximations |14l of those expressions involving Green’s functions (GF) 
that are employed to decide whether there was a reaction or not. This 
is important for 2D systems, because in this case even the radial GF are 
given by complicated integrals that make their numerical evaluation quite 
inefficient. Finally, although the proposed method is designed to provide 
a coarse-grained description, its accuracy can easily be enhanced and the 
method can also be made exact in a natural way. 

Regarding reversible reactions, the algorithm includes possible rebind¬ 
ings after dissociation within the same time step. Simulations show that 
incorporating this effect enhances the accuracy compared to those using 
only a Poisson model. 

The manuscript is structured as follows. First, we give a brief overview 
of those elements of the Smoluchowski theory miEiEnn] needed to explain 
the proposed algorithm. Then, to be specific, we consider a stochastic sim- 
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illation algorithm that has been described before m and that uses exact 
solutions to propagate a molecule pair. Next, we show how the same can be 
achieved much swifter but virtually equally accurate. Finally, we consider a 
few examples and construct, via the proposed simulation method, the PDF 
corresponding to a non-reactive, irreversible and reversible reacting pair in 
3D and 2D, respectively. 


2 Theory 


We consider an isolated pair of molecules with diffusion constants Da,Db 
that diffuse around each other. Alternatively, the pair may be described as 
a point-like particle diffusing with diffusion constant D = Da + Db around 
a static sphere that may be reactive or non-reactive. To be explicit, we 
consider the 3D case, but corresponding considerations apply equally well 
to a molecule pair in 1 and 2D. The PDF that gives the likelihood to find 
the particle located at r = ( x,y,z ), given that it initially was at ro, is the 
GF p(r, t|ro) of the diffusion equation 

^p(r,t\r 0 ) = DV 2 r p{ r,f|r 0 ). (1) 

The GF is subject to the initial 

p(r,t = 0|r 0 ) = <5 (3) (r - r 0 ) (2) 


and BC 

p(r -> oo, t|r 0 ) = 0, (3) 

where r := |r| := \Jx 1 + y 2 + z 2 . The solution to Eqs. (JTJ, ([i]) and (|3j) is 
given by the free-space GF 


p free (r,t|r 0 ) = e 


(r-1-or 


(4) 


Focusing on a radially-symmetric system, the diffusion equation (Eq. 0) 
can be written in terms of the radial coordinate r alone 


d 

—p{ r ,t\r 0 ,t 0 ) = D 


d 2 2d 

Q r 2 r Q r 


p(r,t\r 0 ). 


(5) 


Chemical reactions are incorporated into this formalism by BC at the en¬ 
counter distance. The BC that describe an irreversible and reversible reac¬ 
tion read as follows 
d 

4ira 2 D—p iad (r,t\ r o)\ r=a = K a p rad (r = a, t|r 0 ), (6) 

d 

4ira 2 D —p rev (?’, £|r 0 ) | r=a = K a Prev(r = a, t|r 0 ) - n d [l - S rev (t\r 0 )}, (7) 
or 
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and are referred to as radiation mmm and backreaction BC ia mm m, 
respectively. The radiation BC (Eq. ([6])) describes an irreversibly, partially 
reactive boundary. The limiting cases of a completely reactive (absorbing) 
BC [21] corresponds to fv n -> oo, while a non-reactive pair is described by 
K a = 0 (reflective boundary). The backreaction BC (Eq. 0 ) takes into 
account dissociations also. Clearly, in the limit of a vanishing dissociation 
constant k^ —> 0, the backreaction BC reduces to the radiation BC. The 
survival probability S Tev (t\ro) that appears in the backreaction BC is defined 
by 

/»oo 

Srev(t\r 0 ) = 47T / p re v(r, t\r 0 )r 2 dr. (8) 

J a 

The radial GF corresponding to the different BC are denoted by pf ree , 
Pabsi Prefi Prad and p rev and will play a prominent role in the algorithm we 
will discuss in the following. For later use, we provide the expressions for 
Prad in 3D @1 and p rev in 2D [T5] 


1 


1 



(r - r 0 ) 2 ' 


exp 

ADt 

+ exp 


(r + ro — 2a) 2 
ADt 


Prad(c^o) = 

— kVA irDt exp ( K 2 Dt + (r + ro — 2a) k) erfc ^ K'/Dt + -—^ 

i r°° 2 

p Tev (r,t\to) = — / e~ Dtx 'T(r,x)T(r 0 ,x)xdx, 

^ Jo 

where k := (« a + 47raD)/(47ra 2 D) and the function T(r,x) is defined as 

Jo(rx)f3(x) — Yq (rx)a(x 


(9) 

( 10 ) 


T(r,x ) = 
a(x) = 
/3(x) = 


a 2 (x) + /3 2 (x) 


Kd, 


D 

Kd 


x 2 —— Ji(xa) + 


K n 


2n aD 


xJo(xa), 


x 2 -^mxa) + ^xY 0 (xa). 


( 11 ) 

( 12 ) 

(13) 


2.1 Simulation algorithm 

It is important to emphasize that sampling from any of the radial GF does 
not completely determine the updated 3D positions [9]. Rather, to this 
purpose, one has to use the full solutions describing the angle dependency 
also. The major drawback is that these full solutions are given by quite 
unwieldy integrals, which make their evaluation painfully slow. Although, 
in many situations, the algorithm as a whole is still much more efficient 
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than BD, because large steps can be made, the position update mechanism 
represents a major bottleneck. The analytical representation for the full 
GF describing the diffusion of a point particle around a partially absorbing 
sphere in 3D is known [3] to be 

OO 

PradO, MM = l)P n (cOs(6»)) X 

n=0 

e~ Dtx2 F n+1/2 (r, x)F n+1/2 (r 0 , x)x dx. (14) 

The functions F u are defined by 

p ( x \ — ( 2 ^ + 1 )[J v (rx)Y v (ax) - Y u (rx)J u (ax)\ - 2ax[J u (rx)Yl(ax) - Y u (rx)J' u (ax)\ 
{[(2 h + 1) J v {ax) — 2 axJ' u (ax)] 2 + [(2 h + l)Y u (ax) — 2 axY^ax)] 2 } 1 / 2 

(15) 

Here, 6 denotes the angle between the corresponding relative position vectors 
and one defines h := ha := K a /(47r aD). In the following, we will describe an 
position update method that abandons the use of the full GF altogether. To 
put that update mechanism into context, we focus on a simulation algorithm 
described in Ref. [14] and we briefly summarize its main features. 

We consider a molecule, located at ro = (xo,yo 5 ~o) at time to, close 
to a reactive sphere of radius a. We update the position according to free 
diffusion, i.e. we sample the Gaussian PDF (Eq. Q), which amounts to 
adding the following increments to the Cartesian coordinates 

x = xq + V / 2DAtN(0,1), (16) 

y = y 0 + V2DAfN(0,l), (17) 

z = zq + \ / 2DAtN(0,1), (18) 



where N(0,1) denotes a random number sampled from a Gaussian with 
vanishing mean and variance equal to unity. After this position update, at 
time to + At, the molecule is located at r = (x , y, z). If r = x 2 + y 2 + z 2 < 
a, an encounter took place during the time step At with probability one. 
However, as is well-known, even if r > a, there may have been an encounter. 
To take into account those encounters as well, one may employ radial GF 
M- More precisely, the expression 


p abs (r, Af|r 0 ) 

Pireeif, Af|r 0 ) 


Penc(r, r 0 ,At) 


(19) 


defines the conditional probability P enc that there was an encounter during 
At, given the molecule was propagated according to free diffusion. Thus, to 


5 









test for an encounter, it is sufficient to sample a uniform random number 
£ and to check whether £ < P enc . On the other hand, if there was no 
encounter, nothing else has to be done and the next position update can be 
executed. However, if an encounter was detected, the position r has to get 
corrected (as it corresponds to the free-space GF) and one resamples the 
position according to the expression 


Pref(r, At|r 0 ) - Pabs(r, At|r 0 ), 


( 20 ) 


that involves the full GF (Eq. © with K a = 0 and K a = oo, respectively). 
The rationale to substract the GF with absorbing BC is that p a bs(r> At|ro) 
accounts for exactly those particle trajectories that never encountered the 
boundary during the time step. Finally, the algorithm has to decide whether 
an encounter also led to a reaction. To this end, one employs the conditional 
reaction probability, given that the molecule was propagated according to 


Eq. (20), defined by 


-Preac(r, Hb At) — 1 


Prad(r, At|r 0 ) — Pab S (r, At 17-0 ) 


( 21 ) 


Pref(r, At|r 0 ) -Pabs(r, Af|r 0 ) ' 

This completes a simulation step. As emphasized before, the by far most 


time consuming step is the propagation according to Eq. (20) 


2.2 Update method without full Green’s function 

Now we will show how the time-consuming sampling can be avoided and 
replaced by an iterative scheme that is quite easy to implement. We begin by 
recalling that p re f(r, At|ro) can be numerically approximated by repeatedly 
sampling from Pf ree (r, At|ro) until one obtains r >= a, i.e. all samplings 
resulting in r < a are rejected and redrawn E21 Eg. This sampling scheme 
will result in a renormalized free-space probability density 

Pfree(r, Af|r 0 ) 

/ r > a Pfree(r, At|r 0 )dr' 

Now, as ro —> a or as At —> 0, one has 


Pref(r, Af|r 0 ) ->• 


p fr ee(r, Af|r 0 ) 
/ r > a Pfree(r, At|r 0 )d 3 r' 


( 22 ) 


At first, this seems to indicate that this relation cannot be exploited because, 
first, we would like to keep large time steps and second, in a simulation one 
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faces in general a situation where ro > a and accurately propagating the 
system closer to the boundary requires exactly the measures (either small 
time steps or use of exact full solutions) we seek to avoid. So the central 


question is how we can make use of Eq. (22) when ?’o > a and At is large? 


To this end, we recall that the expressions appearing in the definition of 
the reaction probability (Eq. (|21[)) and in Eq. (20) and that therefore play 
an essential role for propagation and detection of reactions in the algorithm, 
can be written as a convolution relation where the individual factors allow 
for a clear physical interpretation [IT6| 


Pref, rad(Mro) - Pabs(?", t\r 0 ) = 


±D Z 


dT 


f T dr dpabs(r,t- TjQ 

/o T 


Pref, rad(®)T t\Cl) 


<9pabs(£,lizo) 


i=a 


di 


(23) 

£=a 


This decomposition relation motivates to search for an algorithm that is 
capable of faithfully constructing the individual processes represented by 


the rhs of Eq. (23), upon sampling from a free-space density functions only. 


In case of a reflective boundary, sampling according to the renormalized 
free-space density corresponds to the first two factors on the rhs. We now 
detail an iterative bisection method to faithfully model the FP time process 
described by the third factor on the rhs. A similar construction has been 
described in m to determine the last reflection time for a ID problem in a 
different context, but it works equally well in 2 and 3D for our problem at 
hand. 

We reconsider the first part of the previously described algorithm, the 
detection of encounter events via the conditional probability P enc (Eq. ©)■ 
This method is not exact [14] . an error remains due to the ignorance about 
the precise value of the FP time rpp, i.e. the time when r = a for the first 
time. Put differently, the detection of an encounter event via P enc solely 
provides an upper bound, i.e. we only know that to < t fp < to + At. 
However, it turns out that one can determine the FP time and, in addition, 
en passant , the associated full 3D position r rpp with any desired accuracy, 
only upon using P enc and sampling random numbers from the uniform and 
Gaussian densities. 

To see this, we assume that after a position move, the molecule is located 
at r = (x, y, z ) at time to + At. If r < a, there was an encounter. As we 
have discussed before, even if r > a, an encounter may have taken place. 
Hence, we test for encounter by employing P enc as usual. If there was an 
encounter, but r > a, we map the point r outside the sphere to a point 
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inside the sphere by defining 


a — e 

r —y -r, e := r — a. 

a + e 


(24) 


Next, we split the time interval [to, to + At] into [to, to + At/2] and [to + 
At/2, to + At]. Then, we employ the conditional PDF that describes the 
intermediate position ry at time to + At/2, given that the molecule was 
at time to at ro and ended up at time to + At at r. One can show that 
this conditional probability density is again a Gaussian with first moment 
given by (ro + r)/2 and variance cr 2 —> a 2 /2 as follows [13] . Consider two 
independent Gaussian random variables X and Y with variance o 2 and 
vanishing mean. To be definite, we consider the 3D case, but the described 
construction works equally well in 1 and 2D. The quantity we are interested 
in is the conditional PDF p(X = x|Z := X + Y = z). In the context of a 
simulation, the variable Z describes the position update r — ro during a time 
step At, hence X can be interpreted to describe an intermediate position 
update tm — ro during At/ 2. The conditional PDF can be obtained via 


p(x|z) 


P(x, z) 

p( z) 


(25) 


where p(x, z), p( z) refer to the joint PDF of X and Z and the marginal PDF 
of Z, respectively. Using the identities (X 2 ) = (X- Z) = cr 2 , (Z 2 ) = 2cr 2 and 
(X ■ Y) = (Z) = 0, one obtains 


P(x|z) = 


yj (27 t ) 3 <7 3 /2 3 / 2 


exp 


1 (x - z/2) 2 

2 o */2 


(26) 


Hence, it turns out that the conditional PDF assumes the form of a Gaus¬ 
sian with mean z/2 and variance cr 2 /2. This means, in the context of the 
simulation, that an intermediate position can be obtained by a naive BD 
position move (Eq. ©), with adjusted mean and standard deviation). If 
the in this way constructed intermediate point is ym with < a, the FP 
time has to lie between to and to + At/2. But even if ym > a, there may 
have been an encounter in the interval [to,to + At/2]. We test for this as 
usual by invoking p enc . In case there was an encounter, we map the inter¬ 


mediate point inside the sphere (Eq. (24)) and iterate the procedure for the 


interval [to, to + At/2] and for ro and ym- In case there was no encounter 
in that time interval, we have to conclude that the encounter occurred dur¬ 
ing [to + At/2, to + At] instead, and we apply the algorithm to that time 
interval and, correspondingly, to ym and r. Thus, iteratively, we can de¬ 
termine the exact FP time and the exact full 3D position of the molecule 
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when its distance first assumes the value of the encounter radius, only using 
naive BD moves. Next, being located at the encounter radius, we may apply 
the free propagator with rejection to construct the reflective PDF, meaning 
that again we only have to employ sampling from a Gaussian density. Note 
that for this move we have to use the remaining time step, i.e. At — Tpp , 
showing that the first-passage time, that we were initially not interested in, 
plays now an important role. In fact, if one used At, a significant error 
would result. Finally, it is determined whether a reaction took place via 
the reaction probability (Eq. (21)). Note, however, that now one may use 

(r, ro, At). 


-Preac( r - v m , At — Tpp) (compare also with Eq. (23)) instead of P r 


2.3 Reversible reactions 

After the two molecules assumed a bound state, they may dissociate again 
in the case of an reversible reaction. The dissociation probability within a 
time step At is assumed to be given by 


Pdiss(Af) = 1 - exp(—KrfAt). 


(27) 


However, especially when the simulation time step is relatively large, the 
dissociated molecule may rebind again within At. To correct for these rapid 
rebindings, one may use Eq. (21) to take into account the possibility for 
reactions. Thus, the algorithm incorporates dissociations in the following 
way. First, a uniform random number £ is sampled and compared with 
the dissociation probability (Eq. ([27])) . If a dissociation took place, one 
determines the dissociation time tdiss by solving £ = Tdiss(Aiiss)- Then, the 
dissociated molecule is placed at ro = a. The molecule is propagated as 
described previously, where the propagation time step is At — tdiss- Finally, 
one employs P reac for the time span At — tdiss to check whether there was a 
rebinding. 

As we will see in the next section, skipping the test for rebindings, can 
lead to a large error (Fig. [3j right panel). 


3 Simulation results 

We employed the described simulation algorithm to construct numerically 
the GF for the non-reactive and irreversiblly reacting pair in 3D and the 
reversiblly reacting pair in 2D. More precisely, the following simulation set 
up was used: We considered an isolated pair reacting irreversibly A + B —> 
C and reversibly A + B 0 C. Molecule A was held fixed at the origin. 
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Molecule B was placed at an initial position ro = (1.1a, 0, 0). The simulation 
algorithm was run for a time f s im> during which molecule B underwent a 
diffusive motion described by the diffusion constant D. Furthermore, B 
possibly associated with and dissociated from A according to n a and Kd, 
respectively. After each run, we recorded £?’s final postion, unless it was 
bound to A. The corresponding histogram was normalized to account for 
the number of bound states at t s i m . We emphasize that the simulations in 
2D (and At = 0.01a 2 /D) were performed with the small-time expansion 
of the GF p3]. Fig. [I] shows the simulation results for the non-reactive 
and irreversible reacting pair in 3D. The full 3D GF for the non-reactive 
BC (i.e. n a = 0) is shown in Fig. [2j Finally, Fig. [3] shows the results for 
the reversible reacting pair in 2D. All simulation results are compared to 
the corresponding analytical representations (Eqs. @, @, @). We find 
excellent agreement. 
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Figure 1: Numerical construction of the irreversible PDF 47rr 2 p ra( j(r, t|ro) 
for an isolated pair in 3D. The simulation time and the time step are t s i m = 1 
and At = 0.01, respectively. The other parameters are D = a = 1, ro = 1.1. 
The four curves correspond to different values of the association constant 
K a = 0,1,10,100. The markers indicate the height of the histogram and the 
solid lines refer to the exact analytical representation of the irreversible 3D 
GF (Eq. @). 
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Figure 2: Numerical construction of the full GF 2irr 2 sin 6 p Te f(r, 8, t\ro) for 
an isolated, non-reactive pair in 3D. The simulation time is f s ,; m = 1 and the 
other parameters are D = a = 1, ro = 1.1, At = 0.01. The histogram bars 
represent the simulation results and the solid lines refer to the analytical 
representation of the full 3D GF (Eq. (14), with n a = 0). 
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Figure 3: Numerical construction of the PDF 2nrp rev (r,t\ro) for an isolated 
reversible reacting pair in 2D. The simulation time is t s im = 1 and the other 
parameters are D = a = 1, ro = 1.1, n a = 10, At = 0.01. Left panel: 
The four curves describe the analytical results (Eq. (fIo|)) for different values 
of the dissociation constant Kd = 0.01,0.1,0.5,1. The histogram markers 
represent the simulation results. Right panel: Simulation results (Kd = 0.5) 
obtained by using a Poisson dissociation mechanism only (Eq. ©) show 
deviations from the analytical curve (Eq. ©), while results obtained by 
additionally including possible rebinding effects show good agreement. 
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